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ABSTRACT 

This paper develops a novel numerical method for obtaining structures of rapidly 
rotating stars based on a self-consistent field scheme. The solution is obtained iter¬ 
atively. Both rapidly rotating barotropic and baroclinic equilibrium states are calcu¬ 
lated self-consistently using this method. Two types of rotating baroclinic stars are 
investigated by changing the isentropic surfaces inside the star. Solution sequences of 
these are calculated systematically and critical rotation models beyond which no ro¬ 
tating equilibrium state exists are also obtained. All of these rotating baroclinic stars 
satisfy necessarily the Bjerknes-Rosseland rules. Self-consistent solutions of baro¬ 
clinic stars with shellular-type rotation are successfully obtained where the isentropic 
surfaces are oblate and the surface temperature is hotter at the poles than at the equa¬ 
tor if it is assumed that the star is an ideal gas star. These are the first self-consistent 
and systematic solutions of rapidly rotating baroclinic stars with shellular-type rota¬ 
tions. Since they satisfy the stability criterion due to their rapid rotation, these rotating 
baroclinic stars would be dynamically stable. This novel numerical method and the 
solutions of the rapidly rotating baroclinic stars will be useful for investigating stellar 
evolution with rapid rotations. 
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1 INTRODUCTION 

One of the most fascinating challenge in the stellar astro¬ 
physics is understanding the structures of reali stic rotating 
stars and their evolu tion (for example, see IMaeder & Mevnetl 
l200d : [Langerll20 l~3) . According to recent observations, many 
stars display detectable rotations. In fact, massive stars in par¬ 
ticular, suc h as B and O stars , tend to have clearly rapid ro¬ 
tation (e.g. iHunter et al.l[2008l) . Approximately 10 per cent of 
them have the rotatio nal velocities whic h are la rger tha n 30 0 
kms^ 1 teamirez-Agudelo et al.1 120I3I : Dufton et alj l2013h . 
Such rapid rotational velocities can reach the critical veloc¬ 
ity beyond which mass shedding due to rotation occurs and 
no equilibrium state exists oMaederl 120091) . because the typi¬ 
cal values of the critical velocities for 20 Mq star with so¬ 
lar metallicity are ~ 300 kms" 1 - 600 kms -1 during the 
main sequence phase l lEkstrbm et al.l2008l) . The shape of these 
rapidly rotating stars is oblate due to the centrifugal force. 
Therefore, their outer layers cannot be well described by 
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spherical models due to their fast rotation and so non-spherical 
models are required to describe them. The structures of realis¬ 
tic rotating stars are described by using many realistic physi¬ 
cal properties and processes such as viscosity, inhomogeneous 
chemical composition and meridional flows (not considered in 
this paper). In general, the surfaces of equal density (isopyc- 
nic surfaces) and equal pressure (isobaric surfaces) do not co¬ 
incide due to these physical properties and processes, i.e. the 
stars becomes baroclinic. 

Fig. Q] displays two examples of rapidly rotating baro¬ 
clinic stars. The isopycnic surfaces (solid curves) are inclined 
to the isobaric surfaces (dashed curves). The isopycnic sur¬ 
faces are more oblate than the isobaric surfaces in the left 
model, while the isobaric surfaces are more oblate than the 
isopycnic surfaces in the right model. Therefore, the value of 
the density on the isobaric surface varies in the baroclinic star. 
On each isopycnic surface, the density in the left model in¬ 
creases from the pole to the equator, while the density in the 
right model decreases from the pole to the equator. These two 
models are different classes of baroclinic stars explained in 
detail in Sec. 2.1. 
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The baroclinicity is characterized by the angle between 
isopycnic surfaces and isobaric surfaces. It is assumed that the 
system is stationary and axisymmetric and the star does not 
have magnetic field and meridional flow, the Euler equation of 
the rotating star is 


-Vp = -V0 + RQ, 2 e R , (1) 

P 


where p, p, <j> and f l are, respectively, the density, the pressure, 
the gravitational potential and the angular velocity of the star. 
Cylindrical coordinates ( R , ip, z) are utilized in this equation. 
The <p component of the curl of this equation is 

dRfl 2 


1 


(Vp x Vp) = 


dz 


( 2 ) 


If the isopycnic surfaces and the isobaric surfaces coincide, i.e. 
the pressure depends on the density only, the star is barotropic 
(p = p(p)) and the left-hand side of the equation vanishes. 
The rotation of the barotropic star is rigid (fi = flo, constant) 
or cylindrical (fi = because the centrifugal force can 

be derived from a rota ti onal potential (e.g. Eriguchi & Muellen 


Il985l: lHachisull 1986allbl lMaedeJl2009l: lFiiiisawall20l'i) 

r) 7 ? 0 2 

0 = => V^rot. = m 2 {R)e R , 

az 


(3) 


where <f> TO t. is the rotational potential that is an arbitrary func¬ 
tion of R. If the isopycnic surfaces are inclined to the isobaric 
surfaces, the left-hand side of equation 0 does not vanish. 
The angular velocity distributions in the baroclinic stars are no 
longer constant and cylindrical due to the differences between 
isopycnic surfaces and isobaric surfaces. The situation is com¬ 
pletely different from baroclinic star. Therefore, baroclinicity 
is an important component to consider in order to investigate 
such realistic rotating stars and their evolution. 

Nevertheless, only a few works have succeeded in ob¬ 
taining baroclinic equilibrium states with rapid rotation self- 
consistently because it is difficult to calculate and analyse 
baroclinic stars. If it is assumed that the star is barotropic, the 
rotation becomes cylindrical and the centrifugal force can be 
derived from the potential in equation 0. Therefore, the first 
integral of the Euler equation can be obtained analytically as 


/ 


dp 


— -~4> + (plot. + C, 


(4) 


where C is the integral constant. Since the system is barotropic 
(p = p(p)), the integral of the left-hand side of this equation 
can be also calculated. This first integral of the Euler equation 
is used by many previou s works for obtaining rapidly rotat¬ 
ing ba rotropic stars (e.g. lEriguchi & Muelleill 1 9851 : lHachisul 
11986allbi) . In contrast, the situation of the baroclinic star is 
more complicated. If the star is baroclinic, the first integral of 
the Euler equation cannot be obtained analytically because the 
centrifugal force cannot be derived from a rotational potential 
and the distributions of the angular velocity must be calculated 
self-consistently by solving equation 0. 

In order to avoid this difficulty, many previous 
works treated rotating baroclinic sta rs as slowly ro tating 


works treated rotating baroclinic sta rs as slowl y r otating 
stars by using perturbative methods ( Roxburgh et all 1963; 
Roxbur gh & Strittmattedll966 :lFaulkner et al.ll 19681; IClement 

196S; Sackmann^Anand|_ 1970|; Kippenhahn & Thomas 

1970: iMonaghanl 1 197 ll: ISharp et al.lll977l) . Others calculated 


self-consistent rotating equilibrium states, but treated rotating 


stars as pseudo-barotropic, i.e. isobaric surfaces and isopy cnic 
surfaces are not inclined (Jackson 1970; Papaloizou & Whelanl 
1 19731 ; lEriguchi & Muelleil Il99ll ; Ijackson et all 120051) . Al¬ 
though these studies considered the thermal balance of a ro¬ 
tating star in these works, the isothermal surfaces are also iso¬ 
baric and so the system is essentially barotropic. Only a few 
works succeeded in obtaining rapidly rotating baroclinic stars 
self-consistently. 

lUrvu & Eriguchil d 19941) investigated self-consistent ro¬ 
tating baroclinic equilibrium states. They discretized all basic 
equations of fully radiative stars and calculated them simulta- 
neously by using Straight- Forward Newton-Raphson scheme 
dEri smell i & Muellerl 1985). The numerical code was extended 


bv lUrvu & Eriguchil 19951) and more realistic stars with a con¬ 


vectiv e core and radiative envelope were calculated. lRoxburehl 
( 2006 ) obtained self-consistent rotating bar oclin i c stars w ith 
shellular rotation using an iterative method ( Ro xburghll2004) . 
iRoxburghl d2006h did not adopt a particular equation of state 
and did not solve a thermal balance equation. Instead, both 
density and pressure distributions were obtained indepen¬ 
dently and calculated iteratively by fitting the radial profile 
of the angular velocity. [Espinosa Lara & Rieutordl J2007t) suc¬ 
ceeded in obtaining self-consistent rotating baroclinic stars 
with viscos ity and m eridional flows by using a pseudo-spectral 
method fcf. lRieutorcl 1 9871 ). In their numerical code, the phys¬ 
ical variables are projected on to the spherical harmonics and 
solved iterativ ely via a generalized a lgorithm for polytropic 
stellar models dBonazzola et al.lll998l) . More recently, the nu¬ 
merical code was improved and man^_realistic and comple x 
solutions were obtained by [E spin osa Lara & Rieutordl d2013h . 
Very recently. 1 Yasutake et ak ( 20151 ) developed a new formu¬ 
lation and numerical method based on the Lagrangian vari¬ 
ational principle to obtain rotating self-gravitating configura¬ 
tions. They were able to obtain both rotating barotropic and 
baroclinic equilibrium sates consistently by using their new 
method. 


These studies have improved self-consistent rotating 
baroclinic models, but the basic understanding of rotating 
baroclinic stars is still unclear because there are only the few 
solutions of such stars that were developed in the aforemen¬ 
tioned works. Many solution sequences and systematic stud¬ 
ies are required to improve the understanding of rotating baro¬ 
clinic stars and their evolution. Simple and powerful numeri¬ 
cal methods are needed to systematically obtain solution se¬ 
quences. As a first step towards realistic rotating stars and 
their evolution, this paper investigates a versatile numerical 
method for systematically obtaining structures of rapidly rotat¬ 
ing baroclinic stars as well as solutions for baroclinic rotating 
stars. A simple equation of state is adopted in this paper in or¬ 
der to obtain rotating baroclinic stars easily and systematically. 
Some previous self-consistent works adopted complex ther¬ 
mal balance equations, but the rotating baroclinic equilibrium 
state itsel f is determined without any explicit th ermal balance 
equation dRoxburghll2006l ; lYas utake et all2015l ). The mechan¬ 
ical structures themselves are determined independently from 
assumed or additional thermal structures which are totally in¬ 
dependently and freely introduced to mechanical structures in 
Newtonian gravity. It should be noted that the mechanical part 
and the thermal part cannot be separated in general relativistic 
situations because the total energy density is one of the source 
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Figure 1. Seven arbitrary isopycnic surfaces (solid curves) and isobaric surfaces (dashed curves) of baroclinic rotating stars. The outermost curves 
denote the stellar surface. The isopycnic surfaces are more oblate than the isobaric surfaces in the left-hand panel, while the isobaric surfaces are 
more oblate than the isopycnic surfaces in the right-hand panel. These two models are different classes of baroclinic stars. 


of the gravity in the Einstein equations iYoshida et al.ll201~2 ). 
Therefore, in the Newtonian gravity framework, we can ob¬ 
tain rotating baroclinic stars systematically even when we ig¬ 
nore the thermal balance equation. This paper adopts a simple 
equation of state as a baroclinic equation of state, but a real¬ 
istic equation of state could also be easily integrated into the 
method. This new method can calculate rapidly rotating baro¬ 
clinic stars easily. 

The remainder of this paper is organized as follows. Sec¬ 
tion 2 describes the formulation of the rotating baroclinic star 
and numerical scheme of the new method. Section 3.1 veri¬ 
fies the accuracy of the numerical results and compares them 
with results obtained via another numerical scheme. The nu¬ 
merical results and solution sequences are displayed in Section 
3.2. Discussion is given in Sections 4.1 and 4.2. The paper is 
summarized in Section 4.3. 


2 FORMULATION AND NUMERICAL METHOD 

The formulation and details of the numerical method are de¬ 
scribed in this section. The formulation of rotating baroclinic 
stars is summarized in Section 2.1. The remainder of this sec¬ 
tion focuses on explaining the new numerical scheme briefly. 
Both spherical polar coordinates (r, 9, ip) and cylindrical coor¬ 
dinates (R, ip, z ) are utilized in the formulation, but spherical 
polar coordinates are employed in the actual numerical code. 


2.1 A formulation of rotating baroclinic stars 

It is assumed that the system is in a stationary state and the 
configuration is axisymmetric about the rotational axis. It is 


also assumed that the configuration is equatorially symmet¬ 
ric. Both meridional flow and magnetic field are neglected in 
this paper. Although the meridional flow plays importan t role 
for secular transport processes in stellar interiors (e.g. iMathisI 
120131) . the typical time-scale of the meridional flow is much 
longer than the dynamical time-scale of the star. It is almost 
negligible to consider the mechanical structures of the rapidly 
rotating star in equilibrium. Under these assumptions, the ba¬ 
sic equations that describe the rotating baroclinic star are the 
Euler equations 


1 dp 
p dr 


d(j) 

dr 


+ r sin 2 9Q 2 , 


(5) 


1 dp 
~p~d6 


d(j> 

He 


+ r 2 sin 9 cos 9C1 2 , 


the integral form of the Poisson equation 



Pi r ') 

| r — r'\ 


.3 / 

a r 


and a baroclinic equation of state 


P = p(p,X), 


(6) 

(V) 

( 8 ) 


where p, p, </>, f 1, G, and X are, respectively, the density, the 
pressure, the gravitational potential, the angular velocity, the 
gravitational constant and a scalar function that characterizes 
a baroclinicity such as temperature T or specific entropy s. In 
this paper, the following analytic form of a simple equation of 
state is adopted: 


p(p,K) = K(r,9)p 1+1 ' N 
K(r, 9) = K 0 (1 + e 


sin 2 9 i cos 2 9 

1 tf— 
“0 °o 


(9) 
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where e, a o, bo, Ko, m are parameters. K is related to the 
specific entropy of the star. If we assume the ideal gas star 
(p ~ Tp), then K is proportional to the specific entropy 
(K oc exp(s)) and the distribution of K represents the distri¬ 
bution of the specific entropy. If the K at the equator is lager 
than at the pole, the surface temperature is hotter at the equator 
than at the pole. If the K at the equator is smaller than that at 
the pole the surface temperature is cooler at the equator than at 
the pole. Although this equation of state is simple, two types 
of rotating baroclinic stars in Fig.|T]are obtained by this equa¬ 
tion of state. Spherical isentropic (constant entropy) and oblate 
isentropic models were investigated by fixing the values of ag 
and b 0 of this analytic form of K, noting that the isentropic 
surfaces becomes spherical when ao = bo and oblate when 
ao > bo respectively. The spherical isentropic model has the 
spherical isentropic surfaces and the oblate isentropic models 
has the oblate isentropic surfaces. These two types of isen¬ 
tropic models become two types of rotating baroclinic models 
in Fig. Q] The value of e denotes the degree of baroclinicity. 
If we fix e = 0, the system becomes barotropic since K be¬ 
comes constant and p is a function of p (p = p(p)). How¬ 
ever, if we choose e ^ 0, then K is no longer constant and 
the system becomes baroclinic. As the absolute value of e in¬ 
creases, the baroclinicity of the star also increases. As we will 
see in Section 3, this simple equation of state systematically 
provides various rotating baroclinic stars by fixing the parame¬ 
ters in this ways. Since the isentropic surfaces are oblate in the 
oblate isentropic model, the K at the equator is smaller than 
at the pole. Therefore, the surface temperature of the oblate 
isentropic model is hotter at the pole than at the equator. From 
a observational point of view, the surface temperatures of the 
rapidly rotating early type stars are hotter at the p oles than that 
at the equator (e.g. a Leonis observed bylChe et alj|201 ll) by 
the gravitational darkening dvon Zeipelll 19241) . Therefore, the 
oblate isentropic model might be favoured in massive stars, but 
both spherical and oblate isentropic models were calculated to 
obtain various rotating baroclinic stars systematically in this 
paper. 

The ip component of the curl of the Euler equation is used 
in this numerical method instead of the 9 component of the 
Euler equation (equation^, i.e.. 


2 • a a dQ 

r smPcostl—— 
or 


- rsm e ~de 


f dp dp 
\<90 dr 


dp dp\ 
dr d9 ) 


Using cylindrical coordinate, this equation can be rewritten as 

R dj]j_ = ±(dp_dp_dpdp\ 

dz p 2 \dR dz dzdRj ' 1 ’ 


Therefore, we can solve the distribution of Q. by imposing a 
boundary condition on Q at the equatorial surface and inte¬ 
grating these equations directly. In this paper, we impose the 
following boundary condition on at the equatorial surface: 

0 (r,7r/2) = (1 + r 2 /A 2 ) 2 ’ (12) 

where j p, A are constants . This is a j-constant rotation law 
dEriguchi & Muelleil 19851) . 

The physical quantities are transformed to dimensionless 
ones for the actual numerical computations using the central 


density p c and the equatorial radius r e as follows: 

„ r 


t = <t> 

9 ~ 4nGr 2 p c ’ 

■\J 4:7YG pc 


(13) 

(14) 

(15) 

(16) 


k - 


K 

4nGr 2 


l-l/N 


(17) 


where hat'denotes a dimensionless quantity. Note that the di¬ 
mensionless value of the central density is unity in this form 
since 


p c = — = 1. (18) 

pc 

The dimensionless pressure is obtained by the dimensionless 
form of the equation of state as 


P = K 0 


1 + e 



+ 


cos 2 6 \ 



p i+i/n 


(19) 


where e, ao and bo are dimensionless parameters. Thus, the 
central pressure p c can be obtained as 

p c = k 0 pl +1/N . (20) 


The value of ko is determined to ensure a unity dis- 
tan ce from the centre to the equatorial surface of the star 
dTomimura & Eriguchil |2005[) . All equations are also trans¬ 
formed to be dimensionless. For instance, the profile of the 
angular velocity (equation [l2j is transformed as 


H 2 (r,7r/2) 


Jo 


(1 +f 2 /A 2 ) 2 


( 21 ) 


where jo and A are dimensionless parameters. A denotes the 
degree of the differential rotation. 

In summary, the dimensionless forms of equations 0.0 
(lOfind ( 1101 have been defined with the dimensionless equation of 
state (equationll9t for the actual numerical calculation. 


2.2 Numerical scheme 

Next, the numerical scheme of the new method are described 
briefly. The details of the numerical scheme is described in 
Appendix A. The numer ical sc heme ad opts a self-consistent 
field iteration scheme dOstriker & Markl] 19681) . The gravita¬ 
tional field and all physical quantities are calculated iteratively 
in the self-consistent field scheme. In particular, the numeri¬ 
cal schemed of the new method is based on Hachisu’s Self- 
Consistent Field (HSCF) scheme dllachisull 1 986al lhh. which is 
a well-established numerical method for calculating rapidly 
rotating barotropic equilibrium states. However, the HSCF 
scheme cannot calculate rotating baroclinic stars because it 
uses the first integral of the Euler equation (equation 0 . In 
contrast, the numerical scheme of the new method uses the 
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differential forms of the Euler equations (equations l5l II 01 1 and 
can calculate rotating baroclinic stars self-consistently. 

The new method calculate these equations explicitly (Ap- 
pendi x A), while the previous studies for rotating baroclinic 

stars lUrvu & Eriguchill994l.ll995i:lEspinosa Lara & Rieutord 


I2007I . 201 30 adopted implicit calculations (e.g. Newton Raph- 

son method) or spectral-type method to solve these equations 
numerically. Since the implementation of the explicit method 
is much easier than those of the implicit method and spectral- 
type method, the implementation of the new scheme is much 
easier than those of the previous works. We can calculate rotat¬ 
ing baroclinic stars easily and stably by using the new scheme. 

Following the previous works, the axis ratio is 
fixed during iteration cycles in the n umerical scheme 
l lEriguchi &Muelleill 19851: lHachistJfT986allbh . The axis ratio 
q is defined as 


_ t’pol. 


( 22 ) 


where r eq . and r po i. are the equatorial radius and the po¬ 
lar radius of the star, respectively. Extremely rapidly rotat- 
ing stars can be obtained stably by fixin g the axis ratio q 
l lEriguchi & MLiel lerll 1 9851 : lHachisulfl986al) . 


3 NUMERICAL RESULTS 

Numerical results obtained by our novel method are displayed 
in this section. The accuracy of the method is verified in Sec¬ 
tion 3.1.1. In Section 3.1.2, numerical results are compared 
with results obtained via previous method to confirm the nu¬ 
merical results. Two types of rotating baroclinic sequences are 
investigated by fixing the functional forms of the baroclinic 
equation of state with these new solutions displayed and dis¬ 
cussed in Section 3.2. 


3.1 Accuracy verification 

3.1.1 Virial test 


Numerical accuracy and convergence checking are important 
when investigating a new numerical method. To check the nu¬ 
merical accuracy, a relative value of the virial relation was 
computed as follows: 


VC = 


\2T + 3n + W\ 

\w\ 


(23) 


where T, IT and W are, respectively, the rotational energy, the 
total pressure of the star and the gravitational energy with 



(24) 



(25) 


IT = y pdV. (26) 

A numerical domain is defined as 0 < f < 2 in the radial di¬ 
rection and 0 < 8 < 7t/2 in the angular direction. The number 
of the grid points in the r direction within r £ [0,1] and that in 
the 8 direction within 8 £ [0,7r/2] are defined as N r + 1 and 


Ng, respectively. The numerical domains are discretized into 
mesh points with equal intervals A r and A 8 (equations I A11 
and lA2b . Therefore, the total mesh number is 2N r x Ng + 1 
in the computations. Ng = 129 was fixed during the test com¬ 
putations and the value of N r was changed. Both barotropic 
(e = 0) and baroclinic (e = —0.05) solutions were calcu¬ 
lated with q = 0.9, ao = bo = 1-0 and N = 1.5. The 
same model was also calculated using previous codes based 
on HSCF schem e (Fuj isawa et alJ[2012 : lFuiisawa & Eriguchil 
l2013l:lFuiisawjJl2015t) for comparison. 

Fig. [2] displays the convergence of the numerical results 
obtained via each numerical method. The lines denote the ro¬ 
tating baroclinic solutions using the new method (solid line), 
the rotating barotropic solutions using the new method (dashed 
line) and the rotating barotropic solutions using the HSCF 
scheme (dotted line). As seen in Fig. [2] the values of VC be¬ 
come smaller as the grid number in the r direction increases. 
Both the barotropic (dotted line) and baroclinic (dashed lined) 
solutions obtained by the new method converged well. Al¬ 
though the values of the VC obtained by the new method 
are slightly larger than that obtained by HSCF scheme, they 
are sufficiently small and the solutions converged well (see 
lHachisull 1 986al lhl). Therefore, the numerical results by the new 
method are accurate and the numerical accuracy can be consid¬ 
ered verified. In actual numerical calculations, which are tab¬ 
ulated and displayed in this paper, N r = 512 and Ng = 257 
are used. 


3.1.2 Comparison with HSCF scheme 

Next, the solutions obtained by the new numerical method 
were compared with those obtained by HSCF scheme to con¬ 
firm the results. The baroclinic parameter was set as e = 0 
to calculate the rotating barotropic equilibrium states by the 
new method. The same parameters (N = 1.5 and A = 0.9) 
were used and seven models of rotating barotropic equilibrium 
states were calculated via both the new method and HSCF 
scheme. The seven solutions by each code are tabulated in Ta¬ 
ble |T] The left part of the table contains the numerical results 
obtained by the new method and the right one contains those 
obtained by HSCF scheme. The in the table denotes the 
critical rotat ion model beyond which no rotating equili brium 
state exists l lEriguchi & Muelletll 19851 : lHachisJll986ah . The 
isopycnic surfaces of the critical models are displayed in Fig. 

m 

As seen in Table Q] the numerical results of the new 
method (left-hand side) are almost the same as those of the 
HSCF scheme (right-hand side). Although the values of jo are 
slightly different the values of Ko and T/\ W\ are almost iden¬ 
tical. The critical rotation model was also obtained using the 
new method. Fig.[3]shows the isopycnic surfaces of the critical 
rotation models obtained by the new method (left-hand panel) 
and the HSCF scheme (right-hand panel). These isopycnic sur¬ 
faces are almost identical as seen in the figures. The solutions 
of the critical rotation models have cusp structure at the equa¬ 
torial surface (see right-hand panel). The cusp structure is cal¬ 
culated correctly by the new method (left-hand panel). Thus, 
the new method can calculate rapidly rotating barotropic stars 
correctly and accurately. The reliability of the new method is 
thus confirmed by this comparison with HSCF scheme. 
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Figure 2. The values of VC versus the number of the grid points in the r direction ( N r ). Each line denotes rotating barotropic solutions by the new 
method (dashed line), rotating baroclinic solutions by the new method (solid line) and barotropic solutions by the HSCF scheme (dotted line). 


Table 1. Physical quantities obtained by applying the new numerical method (left-hand) and the HSCF scheme (right-hand). The same parameters 
in both codes were used in both codes (N = 1.5, A = 0.9). The solutions with denote the critical rotation models beyond which no rotating 
equilibrium state exists. 


<7 

Ko 

Jo 

T/\W\ 

VC 

New numerical method 

0.900 

2.691E—2 

1.853E—2 

2.232E—2 

7.608E—6 

0.801 

2.368E—2 

3.617E—2 

4.698E—2 

7.712E—6 

0.699 

2.020E—2 

5.255E—2 

7.492E—2 

7.925E—6 

0.600 

1.662E—2 

6.553E—2 

1.048E—1 

8.262E—6 

0.500 

1.272E—2 

7.234E—2 

1.359E—1 

9.087E—6 

0.400 

7.737E—2 

5.616E—2 

1.35 IE— 1 

1.460E—5 

0.395* 

7.126E—2 

5.107E—2 

1.241E—1 

1.658E—5 


3.2 Rapidly rotating baroclinic equilibrium states 

Finally, the rapidly rotating baroclinic equilibrium states were 
calculated using the new method. Two types of rotating baro¬ 
clinic models were developed: a spherical isentropic surfaces 
model and an oblate isentropic surfaces model. The solu¬ 
tion sequences of these two models were calculated system¬ 
atically and the critical rotation solutions of these models 
were obtained via the new method. The Bjerknes-Rosseland 
rulejq which must be satisfied by rotating baroclinic stars (c.f. 
iTassoullf19781) , were checked, and all of the solutions were 
found to them. 


3.2.1 Spherical isentropic surfaces model (do = bo) 

First, the rotating baroclinic stars with spherical isentropic sur¬ 
faces were calculated as one of the simplest baroclinic mod¬ 
els. The parameters of the equation of the state were fixed as 
do = bo = 1.0 to make a spherical isentropic surfaces. Two 
solution sequences (solutions (si-) and (s2-)) with different m 
were calculated by changing the value of q, and ten solutions 
(solutions (sl-a)-(sl-e) and (s2-a)-(s2-e)) and two baroclinic 
critical rotation solutions (solutions (sl-f) and (s2-f)) were ob¬ 


<? 

Ko 

Jo 

T/\W\ 

VC 

HSCF 

0.900 

2.691E—2 

1.835E—2 

2.210E—2 

3.929E—6 

0.801 

2.368E—2 

3.595E—2 

4.668E—2 

4.230E—6 

0.699 

2.021E—2 

5.233E—2 

7.458E—2 

4.646E—6 

0.600 

1.661E—2 

6.543E—2 

1.047E—1 

5.234E—6 

0.500 

1.271E—2 

7.225E—2 

1.357E—1 

6.269E—6 

0.400 

7.702E—3 

5.589E—2 

1.345E—1 

1.030E—5 

0.395* 

7.109E—3 

5.093E—2 

1.238E—1 

1.140E—5 


tained. The physical quantities of the numerical results are tab¬ 
ulated in Table [2] The solutions with in Table [2] denote the 
critical rotation models. The structures of the solutions and 
distributions of the physical quantities are displayed in Fig. [4] 
with solution (s2-a) in the left-hand column, (s2-f) in the cen¬ 
tral column and (sl-f) in the right-hand column. From top to 
bottom, the four rows in Fig. [4] correspond to the isopycnic 
surfaces (solid curves) and isobaric surfaces (dashed curves) 
in the first row, the angular velocity isosurfaces in the second 
row, the K isosurfaces (contours) and the value of K (colour 
maps) in the third row and the K isosurfaces (contours) and 
the value of the dimensionless specific angular momentum j in 
the fourth row. The isopycnic surfaces and isopycnic surfaces 
of solution (s2-e) are also displayed in the left-hand panel of 
Fig-Qto check the baroclinicity easily. Note that Fig.|T|shows 
the arbitrary seven isosurfaces and does not indicate the distri¬ 
butions of density and pressure. 

As seen in the isobaric surfaces (solid curves) and the 
isopycnic surfaces in the left-hand panel of Fig. |T] the isobaric 
surfaces are always more oblate than the isopycnic surfaces. 
This tendency is also found in the first row in Fig. [4] (particu¬ 
larly solution s2-f). This baroclinicity comes from the entropy 
distributions inside the star. The third row shows the K iso- 
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New Method 


HSCF 



A 

R 



R 


Figure 3. The isopycnic surfaces of the critical rotation models (q = 0.395) calculated by the new method (left-hand panel) and the HSCF scheme 
(right-hand panel). The outermost curves denote the stellar surfaces. The difference between the two adjacent contours is 1/10 of the maximum 
value of p. 


Table 2. Physical quantities for the sequences with the spherical isentropic surfaces. The solutions with denote the critical rotation models. 



q 

e 

do 

bo 

m 

K 0 

Jo 

T/\W\ 

n/|w| 

VC 

(s2-a) 

0.900 

0.35 

1.0 

1.0 

2.0 

2.300E—2 

1.390E—2 

1.910E—2 

3.206E—1 

8.447E—6 

(s2-b) 

0.801 

0.35 

1.0 

1.0 

2.0 

1.980E—2 

2.61 IE—2 

3.940E—2 

3.071E—1 

8.817E—6 

(s2-c) 

0.699 

0.35 

1.0 

1.0 

2.0 

1.776E—2 

4.136E—2 

6.579E—2 

2.895E—1 

8.816E—6 

(s2-d) 

0.599 

0.35 

1.0 

1.0 

2.0 

1.472E—2 

5.274E—2 

9.342E—2 

2.71 IE—1 

9.285E—6 

(s2-e) 

0.500 

0.35 

1.0 

1.0 

2.0 

1.133E—2 

5.887E—2 

1.212E—2 

2.526E—1 

1.046E—5 

(s2-f) 

0.408* 

0.35 

1.0 

1.0 

2.0 

7.462E—3 

5.009E—2 

1.238E—2 

2.508E—1 

1.559E—5 

(si-a) 

0.900 

0.35 

1.0 

1.0 

1.0 

2.133E—2 

1.369E—2 

1.967E—2 

3.202E—1 

8.354E—6 

(sl-b) 

0.801 

0.35 

1.0 

1.0 

1.0 

1.899E—2 

2.718E—2 

4.174E—2 

3.055E—1 

8.508E—6 

(sl-c) 

0.699 

0.35 

1.0 

1.0 

1.0 

1.639E—2 

4.009E—2 

6.706E—2 

2.886E—1 

8.757E—6 

(sl-d) 

0.599 

0.35 

1.0 

1.0 

1.0 

1.359E—2 

5.062E—2 

9.452E—2 

2.703E-1 

9.239E—6 

(sl-e) 

0.500 

0.35 

1.0 

1.0 

1.0 

1.047E—2 

5.601E—2 

1.219E—1 

2.520E—1 

1.047E—5 

(sl-f) 

0.406* 

0.35 

1.0 

1.0 

1.0 

6.833E—3 

4.662E—2 

1.225E—1 

2.517E—1 

1.607E—5 


surfaces (contours) and the values of K (colour maps). The K 
isosurfaces are spherical due to the functional form of equation 
J 1 9b . If we assume that the baroclinic ideal gas stars, these so¬ 
lutions have spherical isentropic surfaces because the specific 
entropy s is proportional to K and the K isosurfaces denote 
the isentropic surfaces. The distributions of K in Fig. [4] indi¬ 
cate that the specific entropy increases outward and the equa¬ 
torial specific entropy is higher than the polar specific entropy. 
Therefore, the surface temperature is colder at the poles than 
at the equator. 

The second row shows the angular velocity (Ct) distribu¬ 
tions. The equatorial angular velocity decreases outward be¬ 
cause the profile of the angular velocity on the equatorial plane 
is fixed as equation C3 On the other hand, the angular ve¬ 
locity on the rotational axis increases from the centre to the 
surface ( dCl/dz > 0). Moreover, the values of dCl/dz are 
always positive inside the star. The angular velocity distribu¬ 
tion of a barotropic star is always purely cylindrical (Fig. |3), 


but baroclinicity causes the angular velocity distribution of a 
baroclinic star to not be cylindrical. The distribution is deter¬ 
mined by the inclination angle of the isobaric and the isopy¬ 
cnic surfaces of the rotating star (equation EH). According to 
the Bjerknes-Rosseland rules, the value of dQ/dz is always 
positive if the isobaric surfaces are more oblate than the isopy¬ 
cnic surfaces and the surfac e temperature is colder at the poles 
than at the equator dTassou ill 1978ll . All of these rotating baro¬ 
clinic stars satisfy this rule correctly. The fourth row shows the 
K isosurfaces (contours) and the specific angular momentum 
j = (rsin6*) 2 9 distributions (colour map). As seen in these 
panels, the specific angular momentum on each K isosurface 
(s = constant surface) increases as we move from poles to the 
equator. 

The axis ratios q of both critical models are slightly larger 
than that of the rotating barotropic model (Table Q]. The en¬ 
ergy ratios T/\W\ of the critical rotating baroclinic stars are 
slightly smaller than that of barotropic one. The central and 
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Figure 4. Distributions of the physical quantities. The three columns of figures display the distributions of solution (s2-a) (left-hand panels), solution 
(s2-f) (central panels) and solution (sl-f) (right-hand panels). From top to bottom, the four rows of figures correspond to the isopycnic surfaces 
(solid curves) and isobaric surfaces (dashed curves) (first row), the angular velocity isosurfaces (second row), the K isosurfaces (contours) and the 
value of K (colour maps) (third row) and the K isosurfaces (contours) and the value of the dimensionless specific angular momentum j (fourth 
row). The difference between the two adjacent contours is 1/10 of the difference between the maximum value and the minimum value. 
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A numerical method for baroclinic stars 9 


right-hand panels in Fig. Q] shows the critical rotation model 
with m = 2 (central panels) and m = 1 (right-hand panels). 
Although the isentropic surfaces of m = 1 (sl-f) are more 
concentrated than those of m = 2 (s2-f), the angular veloc¬ 
ity distributions are similar. The parameter m does change the 
distributions of the angular velocity very much. Therefore, the 
shape of the isentropic surfaces is more important for deter¬ 
mining the distribution of the angular velocity. 

Summarizing the above, the rotating baroclinic stars have 
spherical isentropic surfaces and the isobaric surfaces are more 
oblate than the isopycnic surfaces due to the entropy distribu¬ 
tion. Such baroclinicity changes the angular velocity distribu¬ 
tions and the value of dQ/dz is always positive inside the star. 
This con dition itself is one of the Bjerknes-Rosseland rule it¬ 
self dTassoullll978l) . Therefore, all of these rotating baroclinic 
stars satisfy the Bjerknes-Rosseland rules. 


3.2.2 Oblate isentropic surfaces models (do > bo) 

Next, a rotating baroclinic star was calculated with oblate isen¬ 
tropic surfaces. The parameters of the equation of the state 
were fixed as do > bo in order to obtain the oblate isentropic 
surfaces inside the star. A solution sequence with m = 2 was 
calculated by changing the value of q, and five solutions (solu¬ 
tions (o-a)-(o-e)) and a critical rotation solution (solution (o- 
f)) were obtained in the case of the oblate isentropic models. 
The numerical results are tabulated in table [3] The structures 
of the solutions and distributions of the physical quantities are 
displayed in Fig. [5] Solution (o-a) (left-hand column), (o-c) 
(central column) and (o-f) (right-hand column) are displayed 
in Fig. [5] The panels in Fig. [5] are presented similarly to Fig. 
[4] The arbitrary isopycnic surfaces and isopycnic surfaces of 
solution (o-a) are also displayed the right-hand panel of[J]to 
check the baroclinicity easily. 

As seen in the isobaric surfaces (solid curve) and the 
isopycnic surfaces (dashed curves) in the right-hand panel of 
Fig.|T|the isopycnic surfaces are always more oblate than the 
isobaric surfaces. This tendency is also found in the first row 
in Fig. [5] (particularly solution (o-c)). This is the opposite sit¬ 
uation to the spherical isentropic models seen in the previous 
subsection because the entropy distributions near the stellar 
surfaces are different from each other. The third row shows the 
K isosurfaces (contours) and the values of K (colour maps). 
The specific entropy in these models also increases outward, 
but the equatorial specific entropy is lower than the polar spe¬ 
cific entropy. These entropy distributions mean that the sur¬ 
face temperature is hotter at the poles than at the equator. 
This distribution of the specific entropy is opposite to that of 
the spherical isentropic model. The distribution of the angu¬ 
lar velocity is also opposite to that of the spherical isentropic 
model. The second row in Fig. [5] shows the angular velocity 
distributions. The values of dCl/dz are always negative be¬ 
cause of the entropy distributions. This tendency is also in 
contrast to the spherical isentropic model. According to the 
Bjerknes-Rosseland rules, the value of <9f l/dz is always neg¬ 
ative if the isopycnic surfaces are more oblate than the iso¬ 
baric surfaces and the surfac e temp erature is hotter at the poles 
than at the equator iTassoul|[l978h . All of these rotating baro¬ 
clinic stars also satisfy the rules. The shapes of the angular 
velocity isosurfaces are interesting because they are clearly 


shellular-type (Fig. [3- In particular, those of solution (o-a) are 
almost perfectly shellular ones. Although their shapes are de¬ 
formed by the rapid rotations, solutions (o-c) and (o-f) also 
have shellular-type rotation. The fourth row in Fig. [5] shows 
the K isosurfaces (contours) and the values of the specific an¬ 
gular momentum j (colour maps). The specific angular mo¬ 
mentum on each isentropic surface also increases as we move 
from poles to the equator in these models. 

Table[3]shows the solution sequence with the oblate isen¬ 
tropic surfaces. The critical rotation model (solution (o-f)) was 
obtained with the oblate isentropic model. The axis ratio of 
solution (o-f) is larger than that of the barotropic critical ro¬ 
tation and the energy ratio T/\ W\ is smaller than that of the 
barotropic one. The axis ratio is also larger than those of so¬ 
lutions (s2-f) and (sl-f). The right-hand panels in Fig. [5] show 
the critical rotation model (solution (o-f)). The shellular-type 
rotation is realized even if the stellar rotation is as rapid as 
critical rotation. 

Summarizing the above, these rotating baroclinic stars 
have oblate isentropic surfaces, and the isopycnic surfaces are 
more oblate than their isobaric surfaces due to the entropy 
distributions. Such baroclinicity makes angular velocity sim¬ 
ilar to that observed for shellular-type rotation, and the value 
of dfl/dz always negative inside the star. This condition is 
also one of the Bjerknes-Rosseland rules (see e.g. iTassoull 
1 197 8l) . Therefore, all of these rotating baroclinic stars satisfy 
the Bjerknes-Rosseland rules. 


4 DISCUSSION AND SUMMARY 

4.1 Stability of rotating baroclinic stars 

Stability analysis is an important problem to investigate for 
realistic rotating stars. Many types of instabilities are known 
to occur in rotating baroclinic stars, but we restrict ourselves 
to the consideration of instabilities arising from the stellar ro¬ 
tation in this paper. In particular, we focus on the dynamical 
instability due to the rotation that occurs because the time- 
scales of the shear instability arising fr om th e differe ntial ro¬ 
tatio n and the thermal instability (e.g. iGoldreich & Schubert! 
1 19671 ) are much larger than that of the dynamical instability. 
We consider stability criteria derived from linear perturbation 
theory for rotating baroclinic stars. 

The stability criteria of rotating stars are characterized by 
the specific angular momentum ( j = QR 2 ) distributions. Ac¬ 
cording to the Hpiland criterion, a baroclinic star in permanent 
rotation is dynamically stable with respect to axisymmetric 
motions if and only if the two following conditions are satis¬ 
fied: (i) the specific entropy s never decreases outward, and (ii) 
on each surface s = constant (specific entropy constant), the 
angular momentum per unit mass (specific ang ular mom entum 
j) increases from the poles to the equator jTassoulll 19781 ). 

The thermal structures are totally independently and 
freely introduced to the mechanical structures in Newtonian 
gravity, because the thermal structures are obtained by assum¬ 
ing and adding thermal balance equations. If the ideal gas is 
adopted as one of the simplest equation of state, the specific 
entropy is determined and becomes proportional to K. As seen 
earlier in this paper, the specific entropy never increases out¬ 
ward in all of the rotating baroclinic solutions in this paper. 
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Figure 5. As per Fig. [4] The three columns of figures display the distributions of solution (o-c) (left-hand panels), solution (o-c) and solution of a 
critical rotation model (o-f), respectively. 
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Table 3. Physical quantities for the sequence with the oblate isentropic surfaces. The solution with denotes the critical rotation model. 



<? 

e 

CLQ 

bo 

m 

k 0 

io 

T/\W\ 

n/|w| 

VC 

(o-a) 

0.900 

0.45 

1.0 

0.65 

2.0 

1.965E—2 

2.557E—2 

3.404E—2 

3.106E—1 

8.861E—6 

(o-b) 

0.801 

0.45 

1.0 

0.50 

2.0 

1.609E—2 

4.343E—2 

6.365E—2 

2.909E—1 

9.212E—6 

(o-c) 

0.699 

0.45 

1.0 

0.43 

2.0 

1.336E—2 

5.253E—2 

8.812E—2 

2.746E-1 

9.698E—6 

(°-d) 

0.599 

0.45 

1.0 

0.41 

2.0 

1.118E—2 

5.613E—2 

1.076E—1 

2.616E—1 

1.051E—6 

(o-e) 

0.500 

0.45 

1.0 

0.40 

2.0 

8.619E—3 

5.366E—2 

1.202E—1 

2.532E—1 

1.295E—6 

(o-f) 

0.455* 

0.45 

1.0 

0.40 

2.0 

7.324E—3 

4.860E—2 

1.154E-2 

2.563E—1 

1.573E—6 


Therefore, all of the rotating baroclinic stars in this paper sat¬ 
isfy condition (i). The profiles of j on s = constant surfaces 
(K constant surfaces) are displayed in bottom panels in Figs. [4] 
and[5] As shown in these panels, the specific angular momen¬ 
tum on isentropic surface never decreases in both barotropic 
models. All of the rotating baroclinic stars in this paper also 
satisfy condition (ii). Therefore all of the rotating baroclinic 
stars satisfy the Hpiland criterion. They are dynamically sta¬ 
ble configurations against dynamical instability due to rota¬ 
tion. Further analysis of baroclinic rotating stars with realistic 
equations of state would be interest, but thermal balance (en¬ 
ergy equation) must be solved to obtain the baroclinic stars 
with realistic equation of state. It is beyond the scope of the 
present paper. We will construct realistic rotating baroclinic 
models using this new numerical method and systematically 
analyse their stability in future works. 

4.2 Shellular rotation 


Many stellar evolution codes for rotating stars have adopted 
shellular rotation as one of their rotation models (e.g. 


Mevnet & Maedei 1997^ Talon et al. 1997 

; Heser et al. 2000; 

Mathis et al. 2004; Maeder & Meynet 

200aT Potter et al. 

2012allbt Takahashi et al. 20141). The first shellular rotation 


model was formulated and derived bv IZahnl d 1992l~) using a per¬ 
turbative approach. When the turbulence is anisotropic with 
stronger matter transport in the horizontal directions than in 
the vertical, shellular rotation is realized JZahnlll992l) . As a 
result, the horizontal-dependence of the angular velocity dis¬ 
appears and the angular velocity depends on only the distance 
from the centre of the star. This kind of shellular rotation was 
adopted as an initial rotation profile o f core -collapse super¬ 
nova simulations (e.g. lYamada & Satoll 1 994) . However, no¬ 
body had successfully investigated self-consistent configura¬ 
tions of rapidly rotating stars with shellular rotation systemat¬ 
ically because of the difficulty of obtaining rotating baroclinic 
stars. 

Some previous works showed self-consistent 
rotating baroclinic st ars _ with non-cylindrical rota¬ 

ti ons JUrvu & Eriguc M, 119941 : llirvu &Eriguchif [ 19951: 
Esp inosa Lara & Rieutordl 2007 ; lEspinosa Lara & Rieutordl 
20131) , but their rotation distributions were far from shellular. 
On the other hand, [Roxburgh! ( 200^) calculated rotating 
baroclinic stars by fixing the distribution of the shellular 
rotation, but neither calculated solutions syste matic ally nor 
obtained critical rotation models. iKiuchi et al] d2Q10l) investi¬ 
gated multi-layered configurations in differentially rotational 
equilibrium. While they did obtain rapidly rotating solutions 
systematically, they calculated rotating barotropic stars only. 


In contrast, this paper has investigated rapidly rotating 
baroclinic stars systematically as in Section 3.3. The oblate 
isentropic surfaces models were found to have shellular-type 
angular velocity distributions. The equatorial surface tempera¬ 
ture of the rotating baroclinic star is hotter than the polar tem¬ 
perature. These kinds of entropy and temperature distributions 
generate shellular-type rotation. Shellular-type rotation is real¬ 
ized even if the star has very rapid rotation. Solution (o-f) is a 
critical rotation model and its distribution of the angular veloc¬ 
ity is shellular-type (Fig[5]l. These are the first self-consistent 
and systematic solutions of the rapidly baroclinic stars with 
shellular-type rotations. These shellular-type rotating stars are 
dynamically stable as they all satisfy the Hpiland criterion as 
seen in Section 4.1. This kind of shellular-type rotation would 
be realized in many rapidly rotating stars. 


4.3 Summary 

A versatile numerical method for obtaining structures of 
rapidly rotating baroclinic stars was investigated. The novel 
numerical method was based on the self-consistent field 
scheme and the solution was obtained iteratively. Both rotating 
barotropic and baroclinic stars were calculated systematically 
by using the new method. The accuracy of the new method was 
verified by checking the convergence of numerical results. The 
relative values of the virial relation were small and converged 
well. Further, the solutions from the new method and those 
obtained by the HSCF scheme, a well-established numerical 
scheme, were compared, and the solutions were very similar. 

Two types of rotating baroclinic stars — spherical isen¬ 
tropic surfaces models and oblate isentropic surfaces mod¬ 
els — were investigated by fixing functional forms of a sim¬ 
ple baroclinic equation of state. Solution sequences of both 
models were calculated systematically and the critical rota¬ 
tion models were also obtained. All of the rotating baroclinic 
stars satisfy the Bjerknes-Rosseland condition. In the spheri¬ 
cal isentropic surfaces model, the isobaric surfaces are more 
oblate than the isopycnic surfaces and the values of dfl/dz 
are always positive throughout the star. In the oblate isen¬ 
tropic surfaces model, on the other hand, the isopycnic sur¬ 
faces are more oblate than the isobaric surfaces and the values 
of dft/dz are always negative throughout the star. The baro¬ 
clinic star with shellular-type rotation was found to be realized 
when isentropic surfaces are oblate and the equatorial surface 
temperature of the rotating baroclinic star is hotter than the po¬ 
lar temperature if it is assumed that the star is the ideal gas star. 
The shellular-type rotation is realized even if the rotation is as 
rapid as possible. These are the first self-consistent and sys- 
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tematic solutions of the rapidly rotating baroclinic stars with 
shellular-type rotations. 

The Hpiland criterion of the rotating baroclinic stars was 
checked, with all solutions in this paper satisfying the crite¬ 
rion. Therefore, these rotating baroclinic star are stable against 
the dynamical instability due to the rotation. Realistic rotating 
baroclinic stars with a realistic equation of state are needed to 
investigate the evolution of the rapidly realistic rotating stars. 
Further works may extend the numerical method and calculate 
realistic rotating baroclinic stars. 
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APPENDIX A: DETAILS OF NUMERICAL SCHEME 

The numerical scheme of the new method is described in de¬ 
tail. The numeric al scheme is based on a self-consistent fi eld 
iteration scheme fostriker & MartJfl968l : lHachisul 1 1986al lbh. 
The gravitational field and all physical quantities are deter¬ 
mined iteratively. A numerical domain of the computation is 
defined as 0 < r < 2 in the radial direction and 0 < 9 < 
7t/ 2 in the angular direction. The numerical domains are dis¬ 
cretized into mesh points with equal intervals A f and Ad, re- 
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spectively as 

fi = i x Ar, A f=-^~, i = 0,1,2,3, ••• , (Al) 

9j = (j- l)xAJ, = J' = 1.2,3,---.(A2) 

Physical variables are defined at the cross point of two mesh 
lines, i.e., 


Pi ,3 =P(fi,Sj). 


(A3) 


Before the first iteration cycle, a one-dimensional spheri¬ 
cal non-rotating barotropic star is used as the initial guess for 
two-dimensional rotating star. The value of Ko is determined 
to make the radius of the ID spherical star unity. The values 
of Ko and the ID star are used for the iteration. 

We start the iteration cycle after the initial guesses are 
obtained. First, we calculate the gravitational potential <j> from 
the initial guess of p as 


tmax coo 

4> = - ^ 2 «(cos 8) / f 2 f 2 i{r,r)dr 

e=o 
/""V 2 

x / P 2 e(cos9') sinO'dO'p(r ,9'), 

Jo 


(A4) 


where P 2 e is a Lgendre polynomial and f 2 e(r, r') is a function 
defined by 


J r™/r 2e+1 ( f>f') 

f 2 t{r,r) = iy ~ 2 i'i~, 2 t +1 (A5) 

Simpson’s scheme is used to calculate the integrals numeri¬ 
cally. (max is fixed as 2^ max = 40 in all numerical calcula¬ 
tions in this paper. 

Next, we calculate the values of jo and Ko - These values 
are determined in order to make the equatorial radius unity and 
the axis ratio q. First, we calculate the value of Ko by solving 
the dimensionless form of the r component of the Euler equa¬ 
tion (equation^ with the equation of state (equation 1 1 9b on 
the rotational axis. The equations on the rotational axis are re¬ 
spectively discretized as follow: 


the Euler equation (equation[5} and the equation of state (equa- 
tionl 1 9b on the equatorial axis. The equations on the equatorial 
plane are discretized as follows: 


Pi,N g ~ Pi-1,N g _ ( pi,N g + Pi-1,N g ) (<j>i,N e ~ 4>i-l,N e ) 
Ar 2 Af 

, (h+fi- 1) ( pi,N„ + pi-1, Ng) (^i,N e +^i-l,JV s ) /AOS 

I r> r> n ) 


Pi,N e = Ko 


1 +£ I J2 


r 'i,N g 


(pi,r 


\l+-kr 


(A9) 


The value of jo is determined in the same manner to make the 
equatorial radius unity. 

After the new values of Ko and jo have been obtained, 
the distribution of the angular velocity can be calculated us¬ 
ing the dimensionless forms of equations d 1 Ot and d21 b . Four 
neighbouring mesh points [(i — 1, j), ( i, j), (i — 1, j + 1) and 
( i,j + 1)] are used to solve equation d 1 Ob in this numerical 
scheme. The equation is discretized at the centre of the four 
mesh points (left-hand panel of Fie. lAlk The explicit form of 
the discretized equation is 


-2 • a a (^i,j + Dj j+i) — (^i-1,3 + 

r m sin dm cos 8 m -----——--- - - 

2Ar 

- • 2 a (^i.f+1 + ~ ) 

— rm Sin 9m --- 


2A6* 


{.Pi,3 + 1 + Pi-1,j + 1) — (Pi,j + Pi- l,j) 

2A9 


(Pi,3 + 1 +Pi,j) ~ (Pi-l,j+l + Pi-1,j) 
2A r 

(Pi, 3 + 1 +Pt~l,j+l) ~ (Pi,3 + Pi-1, j) 

2A9 


(Pi,3 + 1 + Pi,j) ~ (Pi-l,j + l + Pi-1,3) 


2Ar 


where f m , 9 m and p m are 


(A10) 


fm = - (fi + fi-1) • 


(All) 


Pi, 1 — Pi-1,1 _ (pi, 1 + pi-1,l) (<j>i,l — pi-l,l) 


Ar 


Ar 


(A6) 


— 2 (0j +1 + ■ 


(A12) 


<«) 

A central difference was adopted in this numerical scheme. 
We start the integration at the centre of the star (f = 0). Since 
we have defined the dimensionless form of the central density 
in equations d 1 8b and d20t . the values of pi,i and pi,i are ob¬ 
tained by solving the discretized equations using the values of 
p c and p c . Then, the values of p 2 ,\ and p 2 ,\ are determined 
using the values of pi : i andpi,!. The integration is continued 
until the density pi, i becomes zero {pi,\ = 0), and we obtain 
the value of the polar radius f po i. after the integration. If the 
value of fpoi. is different from the value of q, we change the 
value of Kg and restart t he integration using a shooting-type 
method dPress et al.ll992i ) until we obtain the value of Ko that 
makes the polar radius (f po i.) q- Next, we calculate the value 
of jo by solving the dimensionless form of the r component of 


pm — — (pi-1,j + pi, 3 + f>i—l,j + l + Pi,j + l) ■ (A 13) 

In this scheme, the value of j is calculated by using the 
values of f ii-i j, and Fig. IA1I displays 

schematic images of this scheme. For simplicity, Fig. lAll has 
assumed mesh numbers N r = 3 and Ng = 3. We start the in¬ 
tegration at point 1 in the left-hand panel of Fig. lAll First, the 
boundary condition on the equatorial plane is set. The values 
of Dq, Di 3 , 1 ) 2,3 an d 1 ) 3,3 are fixed during the iterations. At 
point 1 in the left-hand panel of Fig. [All we obtain the value 
of fli i2 from equation j A1 Ob by using the values of Clg and 
Df , 3 because Dq 3 = Do ,2 — Do at point 1. Then, we move 
to point 2 and calculate Cl 2i2 by using the values of Di, 3 , Cl 2 ,3 
and Di, 2 , and solving the equation. 

Finally, we calculate the distributions of p and p from 
equations 0 and £03 by using the distributions of (j> and Cl 2 . 


pi, 1 = Ko < 1 + e 
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Figure Al. Schematic images of the numerical scheme for f2 (left-hand panel) and p, p (right-hand panel). The thick curve denotes the stellar 
surface. (i,j) denotes the mesh point and each number refers to the sequence of the integration. We start the integration of equations at point 1. 


The two mesh points (( i,j ) and (i — 1, j)) are used in this 
calculation. The equation is discretized as 

Pi,j ~ Pi-1,j _ ( Pi, 3 + Pi-1,j) ~ 4>i- 1,J) 

Af 2 Ar 

T^ + ^sin d 2(PiJ + F-i,)(^h+^ 

2 3 2 2 

and the equation of state is given as 


3f&U) 


Pi,j = K 0 < 1 + e 


sin 9 , 


+ 


cos 6j 






A+TV 


(A15) 


We start the integration at point 1 in the right-hand panel of 
Fig. lAll Since the central density and the central pressure have 
been fixed as po = p c = 1 and po = p c by equations H3 and 
( I A1 5t . we calculate pi ,3 and pi ,3 by equation dA14b . Next, we 
move to point 2 and obtain p 2,3 and P 2,3 by using the values 
of pi ,3 and pi t 3 and equations dA14b and dA15l i. We can ob¬ 
tain the values of pip and ptp by using the values of pi-ip 
and pi-ip . If the value of pip or pip is zero, we search the 
location of the stellar surface r s {9), where r s (9) denotes the 
dimensionless radius in each 9 direction. After calculating the 
location of the surface, we move to point 4 and continue the 
integration. After the integration on the rotational axis (point 
6 in the panel), we have obtained new distributions of p and p 
and one iteration cycle is finished. 

In summary, the numerical iteration cycle of this method 
is as follows. 


method in order to make the equatorial radius unity and the 
axis ratio q. 

(v) Calculate the angular velocity 12 by using the distribu¬ 
tion of tj> obtained in (iii). 

(vi) Solve the r-component of the Euler equation with the 
equation of state and obtain the new distributions of density p 
and pressure p by using the distributions of <j> obtained in (iii) 
and 12 obtained in (v). 

(vii) Check the convergence of the system. Compare the 
new and old density distributions. If the relative differences 
between the new distributions and the old distributions are suf¬ 
ficiently small (e.g. smaller than ~ 10 ~ 4 ), the system is con¬ 
sidered to converge well and the iteration cycle finishes as the 
rotating baroclinic star has been obtained in equilibrium. If the 
system does not converge, return to (iii) and start a new itera¬ 
tion cycle until the system converge. 


(i) Set the parameters ( N , q) and fix the functional forms 
of 12 2 on the equatorial plane and the equation of state. 

(ii) Make an initial guess for the two-dimensional rotat¬ 
ing baroclinic star. Calculate a one-dimensional spherical 
barotropic star for initial guesses of the iteration. 

(iii) Calculate the gravitational potential <j> by using the dis¬ 
tribution of p. 

(iv) Calculate the values of Kq and jo via a shooting-type 
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